%
% Solve Shape Optimization Problem using moving mesh finite element method
%  the elasticity problem is solved using linear fe on triangle
% and the phase field is calculated on the node using linear fe
%

%% 1. physical parameters
e = 1.0; %Young's modulus
nu = 0.3; %Poisson's ratio


%% 2. the computational mesh
lx = 2.0; %Length of x
ly = 1.0; %Length of y
nelx = 100; %Discritization number of x
nely = 50; %Discritization number of y
h = lx / nelx; %Element size
[V,T] = quad_tri(0,lx,0,ly,nelx,nely);
[E,ET,TE,VT,ring1,NodeFlag] = build_geom_info(V,T);

%% 3. computational parameters
lv = [0 -1];  % Local load vector, must consider the element 
lp = [1 1];   % originally at point (nelx+1)*(nely/2+1); 
              % Load point nodes number (center of the right side)
kappa = 1e-4; % Coefficient of diffusion term (The most important parameter!
              % In general, it should be set in the order of dx^2. But, in my method, 
%smaller value than the basis works well.)  

PDECycle = 20; %Number of the time updating in Allen-Cahn equation per one iteration.
%In some sensitive problems (such as compliant mechanism problem or vibration problem),
%it should be small.
timeStep = 0.5 * h^2/(4*kappa);
%It should be desided satisfying CFL condition.
tol = 1e-3; %Tolerance for objective function increasing
cycle = 20; %Number of iteration


%% INITIALIZE
[V,phi] = init_phi_circle7(V,T,VT,ring1,NodeFlag);  %Set initial value of phase field function phi.


%% begin the main cycling
for l = 1:cycle 

    % this is not good 
%    phi =  max(0,min(1, phi)); %Clip phi the range between 0 to 1. This is a charm against numerical errors.

%              
%     colormap(gray);
%     imagesc(flipud(-phi'),[-1.0 0]);
%     axis equal; axis tight; axis off;pause(1);%Above 3 lines are plotting phi. 
%       
    
    % 1. solve the state equation with finite element method
    u = FEM_u(V,T,e,nu,lv,lp,phi); %Calcurate displacement u by FEM.
  
    
    % 2. calculation of the compliance and its sensitivity,
    % need to be modified to accomodate the nonuniform mesh environments
    [com, vol, objf, mphi] = forcing(V,T,u,phi,e,nu);
    
    disp([' It.: ' sprintf('%4i',l) ' Obj.: ' sprintf('%10.4f',objf) ' Vol.: ' sprintf('%6.3f',sum(vol))]);
   
    % Solve Allen-Cahn equation with several time steps
    phi = FEM_phi(V,T,kappa,timeStep,mphi);
        
    % 4. adaptation of the mesh    
%    [V,phi,u] = update_mesh(V,T,VT,ring1,NodeFlag,phi,u);
    
end